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Abstract 



We propose new algorithms for numerical integration of the equations 
of motion for classical spin systems with fixed spatial site positions. The 
algorithms are derived on the basis of a mid-point scheme in conjunction 
with the multiple time staging propagation. Contrary to existing predictor- 
corrector and decomposition approaches, the algorithms introduced preserve 
all the integrals of motion inherent in the basic equations. As is demonstrated 
for a lattice ferromagnet model, the present approach appears to be more 
efficient even over the recently developed decomposition method. 
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The prediction of collective phenomena in magnetic materials was the subject of many 
investigations in theory and computer experiment These investigations dealt mainly 

with static properties, in particular, with phase transitions, critical behavior and scaling. 
The theoretical description of dynamic properties presents a much more difficult problem 
which, at the present time, cannot be solved quantitatively even for the simplest lattice 
systems such as the Ising, XY, and Heisenberg models. Until now, the method of molecular 
dynamics (MD) can be considered as the main tool for quantitative studies of dynamic 
critical behavior, dynamic scaling, and processes of spin relaxation. 

The construction of stable and efficient MD algorithms for systems with spin (orienta- 
tional) degrees of freedom remains a current problem. The traditional numerical methods |7j 
on integrating differential equations are unsuitable because of their high instability on MD 
scales of time. For this reason, standard predictor-corrector schemes were utilized in recent 
years for MD simulations of classical XY and Heisenberg models in d = 2 and of Heisenberg 
ferro- and antiferromagnets in d — 3 dimensions @-l2|. However, the requirement on total 



energy conservation restricts these schemes to be used with very little step sizes only. 

Quite recently, new spin dynamics algorithms have been devised [P^-TT5|]. They are 



based on the Suzuki- Trotter (ST) decompositions of exponential operators and unlike usual 
methods preserve the total energy to within machine accuracy. It was shown that these 
algorithms allow much larger time steps than the predictor-corrector schemes and thus may 
lead to a substantial speedup of MD simulations. However, the decomposition integrators 
destroy the conservation of magnetization vector. Moreover, they are applicable, in fact, 
only to spin systems when the decomposition into noninteracting sublattices is possible and 
cannot be used for models with arbitrary lattice structures. 

In this report we solve the problems just mentioned and derive algorithms which preserve 
all the conservation laws imposed by the equations of motion. The algorithms are tested 
and compared with previous approaches. 

Let us consider a collection of N spins represented by continuous three-component vectors 
Si = (sf,sf,Si) with the fixed length |sj| = 1 for each site i. A typical model Hamiltonian 



for such a system can be cast in the form 

N N 

H = - £ Jistfq + s\s) + \s*s*) - C^) 2 (1) 

i<j i=l 

where Jij is the exchange integral for a pair of spins, A is the exchange anisotropy 
parameter, and C denotes the strength of single-site field anisotropy. At C = 0, Eq. (1) 
represents the isotropic (A = 1) or anisotropic (A ^ 1) Heisenberg ferro- or the corresponding 
antiferro-magnet for J > and J < 0, respectively. For A = C = 0, Eq. (1) reduces to 
the XY model. We do not restrict ourselves to lattice systems with the nearest-neighbor 
interaction, and the results presented below can be used for continuum models with arbitrary 
spatial spin distributions as well. Therefore, we indicate explicitly by the subscripts i,j that 
the exchange integral depends on spatial positions (which are fixed but not necessarily 
periodic) of spin sites. 

The dynamic properties of the system can be obtained from MD simulations by numerical 
integrating the following equations of motion [|S|-|T3|j: 

Si = —s i = ^^-Xs i (t) = n i (t)Xs i (t), (2) 

where fij = — ^(Yljfj^i) Jij( s ji ^ s j) + 2(7(0,0, sf)) denotes the local Larmor frequency. 
Since the effect of collective thermal excitations (e.g., phonons) is not described by the Hamil- 
tonian (1), Monte-Carlo (MC) simulations must be employed additionally |l3j to generate 
equilibrium configurations as initial conditions to Eq. (2). This procedure is justified by 
the fact that in magnetic systems the characteristic time intervals corresponding to varying 
spin variables are much shorter than typical time scales of the thermal excitations. 

In view of the symmetry = Jji, it follows from Eqs. (1) and (2) that the total energy 
E = H is an integral of motion, i.e. dE/dt = 0. The magnetization M = J2i s i is a ls° 
conserved during the spin evolution of the isotropic Heisenberg model. For the anisotropic 
case (A 7^ 1 and/or (7^0) only the component M z of M will unchange in time. In addition, 
the structure of Eq. (2) imposes also the conservation of individual spin lengths. Existing 
MD algorithms do not fulfill these conservation laws simultaneously. Thus, in order to 
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reproduce the dynamical behavior properly it is required that the deviations of conservative 
quantities from their exact values to remain within an acceptable level of precision. This 
leads to obvious limitations on the size of time steps which can be used in MD simulations. 
It would, therefore, be very desirable to derive algorithms which conserve all the integrals 
of motion exactly or, at least, within machine accuracy. 

The basic idea of our approach consists in the following. Suppose that an initial spin 
configuration {sj(t)} has been specified and we would like to obtain values of Sj at time t + r 
within 0{t z ) truncation terms, where r denotes the step size. This can be realized using 
a mid-point scheme, Sj(i + t) — Si(t) + Si(t + t/2)t + 0(r 3 ). The time derivative can be 
determined applying the usual interpolation formula Sj(t + r/2) = |[sj(i) + Sj(£ + r)] + 0(t 2 ). 
Such a formula, however, does not maintain the unit norm of Sj(i + r) and thus needs in 
modifications. Since Sj depends on both the local frequency fij and spin orientation Sj, it 
is more natural to apply the interpolation with respect to these two dynamical variables 
separately rather than to the function Sj as a whole. In doing so we obtain Sj(i + r/2) = 
i[(f2j(t) + Qi(t + t)) x (si(t) +Si(t + r))] + 0(t 2 ), resulting in the implicit spin propagation 

s! n+1) (t + r)=B i (t) + ^[ n<"> X (s, (t) + (t + r))] 

+ 0(r% (3) 

where Q, = |[0;(t) + fl^\t + r)]. As far as the mid-step frequency tl^ depends itself 
on (generally speaking) all spin orientations at time t + r, Eq. (3) constitutes a set of iV 
quadratic equations for {Sj(t + r)} which can be solved iteratively (n = 0, 1, 2, . . .) putting 
s\°\t + r) = Si(t) + fli(t)xsi(t)T as the initial guess (note that iterative solutions are inherent 
in predictor-corrector and decomposition schemes too). Already one iteration is enough to 
reach the required 0(r 3 ) accuracy. The necessity of performing further several updates of 
Eq. (3) will be motivated latter. 

It can be shown readily that for the isotropic model the magnetization is conserved 
exactly during the spin dynamics propagation given by Eq. (3). Indeed, summation of Eq. 
(3) over the spin numbers and taking into account the explicit expression for fl^ yields 
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M ("+ 1 )(t + r) = M(t) + AMW, where AM (n) = ^ £ ¥j Jij^ii*) + sj n) (t + r)] X [sj n) (f) + 
sj n ^(t + r)]. The term AM'"' is canceled because of the invariance of the double sum with 
respect to the transformation i <-> j, and of the obvious equality axb + bxa = which 
is valid for arbitrary vectors a and b. Thus, M(£ + r) = M(i) within each iteration. The 
proof of the conservation M z (t + r) = M z (t) for the anisotropic case is similar. 

Another important feature of the mid-point propagation is that it conserves the total 
energy within machine accuracy 0(e), where e denotes the iterative precision, i.e., \s\ n+1 \t+ 
t) — s\ n \t + r)| < e. To show this, let us perform a scalar multiplication of Eq. (3) with 
the vector f2j = §[f2j(£) + toi(t + r)]. Then using the equality [axb]-a = one obtains 

Si (t + r).[fii(*) + + r)] = Si(*) + + r)], 

where 0(e) terms have been neglected. Summing up the last relation leads to E(t + r) = 
| E l Si(t + T)'ili(t + t) = | E 4 Si^-fi^t) + AE, where AE = | Ejs^i)-^^ + r) - s,(t + 
r)«f2j(t)]. The term AE is canceled again because of the linear dependency of fli on spin 
components, and of the symmetry = Jji, so that E(t + r) = E(t) + 0(e). The uncertainty 
e can be reduced to a negligibly small value at a given r by adjusting the number I > 1 of 
iterations for Eq. (3). The rapid convergence e — > +0 is guaranteed by the power dependence 
e ~ (9(V +2 ) and by the smallness of r. Of course, the iterative solutions require additional 
computational efforts, but they are compensated completely by using larger time steps. For 
instance, spending the same amount of computer time, we could try to reduce the energy 
deviations within only one iteration by decreasing the time step to r/l. This way, however, 
is very inefficient because then the deviations will behave like C((r//) 3 ) and, thus, decrease 
with increasing / much more slower than the power dependence 0(r l+2 ), in other words 
r l+2 < (r//) 3 . 

An additional surprising property of the mid-point integration is the conservation of spin 
lengths. Implicit evaluations given by Eq. (3) achieve this conservation in iterative sense, 
i.e., \Si(t + t)\ = 1 + 0(e). In order to maintain spin lengths exactly, the iteration process 
should be reconstructed. Considering the quantity tt^ as a parameter, Eq. (3) can be 
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solved analytically, 



1 



XSi(t)r 



+^(20S" ) (Or ) .s { ,(t))-(^ ) ) 2 s J (t)));, 



(4) 



and further iterated because of the explicit dependence of Qi on spin orientations {sj(t + r)}. 
Obviously, such a modified iterative scheme will conserve the magnetization and total energy 
(to within machine accuracy) like the usual spin evaluation (3) (since Eq. (4) was obtained 
from Eq. (3) by the exact transformation). Moreover, as can be verified easily, Eq. (4) 
presents the unitary propagation Sj(i + r) = Dj(t, r)Sj(t), where Dj(t, r) is an orthonormal 
matrix which rotates the vector Sj(i) on angle ip = arcsin(f2jr/(l + r 2 Qf /4)) around axis f2j. 
Therefore, the modified scheme will maintain the unit norm of spin lengths perfectly, i.e. 
|sj(t + r)| = \si(t)\ = 1, for each iteration. This may lead to more efficient computations, 
despite a somewhat complicated structure of the RHS of Eq. (4) with respect to that of Eq. 



The convergence of Eq. (4) can be improved significantly by using an advanced iterative 
method. Namely, recalculating a current (i = 1, 2, . . . , N) value of s- n+1 ^(t+r) within the nth 
iteration, it is necessary to take into account the already obtained quantities sjj. n it + r) for 
k — 1, 2, . . . , i — 1 when forming the RHS of Eq. (4). The advanced method is preferable for 
lattice systems with the nearest-neighbor convention, when the computer time per iteration 
required for the recalculations of spin values (according to Eq. (4)) dominates over the time 
needed to update the mid-step Larmor frequencies. Further optimizations are also possible 
for each specific model. This completes our mid-point spin dynamics (MPSD) algorithm of 
the second-order. 

A way to construct higher-order versions of the MPSD integrator lies in employing a 
multiple staging technique used earlier |13[ in the framework of the ST approach. We have 
realized that this technique is applicable for our approach too and obtained the following 
result 



(3). 
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*(* + r) = II D,(t,epT) Si (t) + 0(r m+1 ), (5) 
P =i 

where the coefficients £ p are chosen at a given number P in such a way to provide the highest 
possible value for m. The desired fourth-order (m = 4) algorithm (MPSD4) can be directly 
derived from Eq. (5) using P = 5 and the coefficients £i = £ 2 = £4 = £5 = £ = 1/(4 — 4 1 / 3 ), 
and £3 = 1 — 4£. That is very interesting, these coefficients coincide with those obtained 



within the ST decomposition of exponential operators [|i~6[|. 

Clearly, the fourth-order version is energy- and magnetization-preserving and conserves 
spin lengths (since the conservations are achieved at each stage p). The solutions generated 
by the MPSD/MPSD4 algorithms are also time reversible (because past and future values 
of Si and f2j enter symmetrically into the interpolated function Sj, and because £ p appear 
symmetrically in Eq. (5)). The reproduction of the last feature is particularly important 
as well since the numerical stability of an algorithm is directly connected with its time 



reversibility [17 



In our MD simulations we considered a simple cubic lattice in d = 3 with N = 1000 
sites imposing periodic boundary conditions to each direction. The strongly anisotropic 
case C = J with J > 0, = J5ij and A = 1 was chosen to describe spin interactions 
(Eq. (1)), where 5^ = 1 for the nearest neighbor j of site i, and 5y = otherwise. All 
test runs were started from an identical well equilibrated configuration prepared by us with 
the help of MC simulations at a temperature T = 0.8T C , where T c = 1.442929 J/k^ is the 
critical temperature ]TJ| of the isotropic model (C = 0). The simulations were performed 
on the Origin 2000 workstation at the Linz University. The equations of motion were 
integrated using the Adams-Bashforth-Moulton (ABM) predictor-corrector integrator @ at 
r* = rJ/h = 0.01, the ST decomposition schemes 0] of the second (STD) and forth (STD4) 
orders at r* = 0.04 and 0.2, respectively, as well as using our MPSD and MPSD4 algorithms 
(Eqs. (4) and (5)) at r* = 0.04, 0.1, 0.2, and 0.4. 

Examples on the total energy E* = E/J and magnetization M z conservations are shown 
in Fig. 1. The huge energy drift (see dashed curve in Fig. la) indicates clearly that the ABM 



algorithm is unsuitable for long-duration observations even at the smallest time step. This 
is explained by the irreversibility of the ABM integrator and the fact that it destroys the 
unit norm of spin lengths. At the same time, the STD/STD4 algorithms allow much larger 



step sizes, that is in the self-consistency with a conclusion of Ref. ||13|| . Three iterations 
were sufficient for the STD algorithm to obtain a level of energy conservation presented in 
Fig. la. Correspondingly 5 and 6 iterations were required for the STD and STD4 algorithms 
to conserve the total energy within machine accuracy (e ~ 1CT 9 in our program code). The 
STD/STD4 integrators, however, do not conserve the magnetization (see Fig. lb) which 
fluctuates quite visibly especially in the STD case. These fluctuations are caused by the 
0(t 3 ) (or (9(r 5 )) truncation errors and will increase drastically with further increasing r. 

The pattern is different for the MPSD/MPSD4 algorithms because they can conserve 
both the total energy and magnetization at, in principle, arbitrary time steps. We have 
determined the following values for the number / of iterations needed for obtaining the con- 
servation to within machine accuracy: I = 5,8, 11, and 18 for the MPSD as well as I — 4, 6, 7, 
and 11 for the MPSD4 integrators, corresponding to the time steps r* = 0.04,0.1,0.2, and 
0.4, respectively (note that the advanced method was used to iterate Eq. (4)). The MPSD 
(MPSD4) algorithm required at r* = 0.04 (r* = 0.2) approximately the same computer 
time as that of the STD (STD4) scheme for the integration over the fixed t* = tJ/H = 1000 
interval. However, the MPSD/MPSD4 algorithms can be used with larger time steps 
(r* > 0.04/t* > 0.2) in view of their energy- and magnetization-preserving properties. 
Taking into account that the number / increases with increasing r slower than linearly, these 
algorithms will lead to an improvement efficiency of the computations. Of course, we can- 
not apply too large step sizes (r* ~ 1), because then the microscopic solutions will deviate 
considerably from exact trajectories. The final decision on using the biggest possible values 
of step sizes for the MPSD/MPSD4 algorithms can be done in each specific case by direct 
MD measurements of macroscopic observable quantities. 
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Our measurements were performed for the dynamic structure factor 

S(k,u) = — / (^s, T (0)-sJ(t)e ik -( r -^))e^dt, 

which allows one to extract the spectrum of collective transverse spin excitations, where 
denote the lattice vectors and is a perpendicular to M component of Sj. The function 
S(k,u) is plotted in Fig. 2a in the undimensional presentation S*(k,cu) = JS(k,cu)/h with 
uj* = cuh/J, k* = k/k min and fc min = 2ir/N 1/3 = n/5. It was obtained within the MPSD 
integration at r* = 0.04, and the averaging ( ) was taken over the time t* = 1000 for each of 
1000 runs (with independent initial MC configurations). In order to investigate the influence 
of increasing step size on the spectrum, the calculations were repeated with the STD/STD4 
and MPSD/MPSD4 algorithms at various step sizes. The frequency uj* iax and the relaxation 
time r* or of the transverse waves, calculated numerically at k* — 1 from the peak position 
and its half-width, respectively, are shown in Fig. 2b and 2c as the functions of the time 
step. As can be seen clearly, the MPSD/MPSD4 algorithms are less sensitive to increasing r 
than the STD/STD4 integrators. For instance, the level of deviations in 0J^ ax obtained with 
the STD (STD4) schemes at r* = 0.04 (r* = 0.2) can be observed with the MPSD (MPSD4) 
algorithms at step sizes which are approximately in factor 1.5-2.0 larger, namely, at r* ~ 0.07 
(t* ~ 0.3). This gain in time is explained by additional cancellations of truncation errors 
due to the conservation of all the integrals of motion within our approach. For the same level 
of accuracy, the fourth-order versions MPSD4/STD4 allow step sizes which are nearly in 4-5 
times larger than those of the second-order algorithms MPSD/STD. This compensates to 
some extent the additional computational efforts needed to evaluate high-order expressions. 
However, if very high precision is required, the fourth-order schemes become more efficient, 
because then the truncation errors decrease more rapidly with decreasing the step size. 

In the conclusion we point out that alternative algorithms for classical spin dynamics 
simulations have been proposed. Their advantages over the predictor-corrector schemes are: 
(i) time-reversibility, (ii) exact conservation of spin lengths, and (iii) allowance of much larger 
time steps. The advantages over the decomposition integrators consists in (i) magnetization 
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conservation, (ii) allowance of larger step sizes, and (iii) applicability to systems with an 
arbitrary lattice structure and to continuum spin models. Moreover, the possibility of the 
new algorithms to preserve all the conservation laws should be considered as the chief feature 
which distinguishes them from all existing MD integrators. There are no other algorithms 
of such a kind known for any system of interacting particles. This fact may play a role in 
the methodology of MD and stimulate further investigations on constructing conservation 
laws preserving algorithms for other systems. 

This work was supported in part by the Fonds zur Forderung der wissenschaftlichen 
Forschung under Project No. P12422-TPH. 
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Figure captions 

FIG. 1. The total energy E*/N (subset (a)) and magnetization M z /N (subset (b)) per 
spin as functions of the time length t* = tJ/k of the simulations carried out for a lattice 
ferromagnet model using the predictor-corrector (dashed curve, marked as ABM), decom- 
position (solid curves, STD/STD4), and mid-point (bold solid lines, MPSD) approaches. 

FIG. 2. (a) The reduced transverse dynamic structure factor S(k,u) of a lattice ferro- 
magnet model; (b) and (c): The frequency peak position (b) of S(k*,uj) and the relaxation 
time (c) of transverse waves as functions of the time step r*, obtained at k* — 1 within 
different integration schemes. Numerical errors are less than the size of symbols used in the 
subsets (b) and (c). 
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